Introduction

The goal of study 2 was to…

Methods

library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)
source("functions.R")

df <- read.csv("../data/study2.csv") |>
  mutate(
    Date = lubridate::dmy(Date),
    Participant = fct_reorder(Participant, Date),
    Screen_Refresh = as.character(Screen_Refresh),
    Illusion_Side = as.factor(Illusion_Side),
    # Illusion_Effect = fct_relevel(as.factor(Illusion_Effect), "Incongruent", "Null", "Congruent"),
    Illusion_Effect = fct_relevel(as.factor(Illusion_Effect), "Incongruent", "Congruent"),
    Block = as.factor(Block),
    Education = fct_relevel(Education, "Doctorate", "Master", "Bachelor", "High School", "Other")
  )

# Fix precision
for(ill in unique(df$Illusion_Type)) {
  data <- df[df$Illusion_Type == ill, ] 
  i <- 10
  while (length(sort(unique(signif(data$Illusion_Difference, i)))) != 8) {
    i <- i - 1
  }
  df[df$Illusion_Type == ill, "Illusion_Difference"] <- signif(df[df$Illusion_Type == ill, "Illusion_Difference"], i)
  if (i != 10) {
    message(ill, ": Illusion_Difference values != 8. Rounded to ", i, ".")
  }
}

# Transformation
df <- df |> 
  mutate(    
    Illusion_Difference_log = log(1 + Illusion_Difference),
    Illusion_Difference_sqrt = sqrt(Illusion_Difference),
    Illusion_Difference_cbrt = Illusion_Difference**(1 / 3),
    Illusion_Strength_log = sign(Illusion_Strength) * log(1 + abs(Illusion_Strength)),
    Illusion_Strength_sqrt = sign(Illusion_Strength) * sqrt(abs(Illusion_Strength)),
    Illusion_Strength_cbrt = sign(Illusion_Strength) * (abs(Illusion_Strength)**(1 / 3))
    )

Exclusions

# plot(estimate_density(dfraw[dfraw$Participant == "60684f29dbfe1bb2059e5e27_rkqoy", "RT"]))

# Dear participant, thank you for participating in our studfy. Unfortunately, we didn't receive your data :( did something happen? Some technical issue? We would like to kindly ask you to return your participation so that we can open up more slots. Thank you in advance, and apologies for the inconvenience! 

# Dear participant, thank you for participating in our study. Unfortunately, our system detected multiple issues in your data (such as implausibly short responses, and random-like pattern of answers), which makes it unusable. We understand that you might have been in a hurry or had some other issues; we hope to open-up more slots in the future would you be interested to participate again.

# Dear participant, thank you for participating in our study. Unfortunately, our system detected multiple issues in your data (such as implausibly short responses, and random-like pattern of answers - in particular in the 2nd part of the study), which makes it unusable. We understand that you might have been in a hurry or had some other issues, and so we kindly ask you to return your participation; we hope to open-up more slots in the future would you be interested to participate again.

# Just received the results: in your case, the three most prominent issues were your response pattern that was equivalent to random (in the sequence of keystrokes) that led to a chance-level performance (that was also significantly different from the rest of the population). Moreover, your reaction time distribution was also very different from the norm, with a vast majority of implausibly short responses (i.e., that are faster than the time it takes the brain to process a visual input). This issue was even more prominent in the second block (after the break), which typically happens when participants are in a rush to finish. Finally, your overall completion time was also significantly below the average. Again, we apologize, we understand that your time is valuable, but unfortunately we run too on limited funds and can hardly spare more on unusable data. Sorry for the inconvenience, we will try to open-up more slots soon.

outliers <- c(
  # Error rate of 48.8% Very short RT
  # Prolific Status: REJECTED (06/08)
  "5e66ed53de7896486f9a71f8_p8ckk",
  # 2nd block of responses very fast
  # Prolific Status: REJECTED (15/08)
  "61443e47c248a375c899a0cf_9qz3u"
)

partial_outliers <- c(
  # 2nd block a bit bad
  "5c43cd414fe4f800016e4983_0zp37", 
  # Entire 2nd block bad
  "61564e378e974bdbe42763a2_hhufm")

We removed 2 participants upon inspection of the average error rage (when close to 50%, suggesting random answers) and/or when the reaction time distribution was implausibly fast.

For each block, we computed the error rate and, if more than 50%, we discarded the whole block (as it likely indicates that instructions got mixed up, for instance participants were selecting the smaller instead of the bigger circle).

Error Rate

dfsub <- df |>
  group_by(Participant) |>
  summarize(
    # n = n(),
    Error = sum(Error) / n(),
    RT_Mean = mean(RT),
    RT_SD = sd(RT),
    IPIP6_SD = mean(IPIP6_SD),
    PID5_SD = mean(PID5_SD),
  ) |>
  ungroup() |>
  arrange(desc(Error))

data.frame(Participant = c("Total"), t(sapply(dfsub[2:ncol(dfsub)], mean, na.rm=TRUE))) |> 
  rbind(dfsub) |> 
  knitr::kable() |>
  kableExtra::row_spec(1, italic = TRUE) |> 
  kableExtra::row_spec(which(dfsub$Participant %in% outliers) + 1, background = "#EF9A9A") |> 
  kableExtra::row_spec(which(dfsub$Participant %in% partial_outliers) + 1, background = "#FFECB3")
Participant Error RT_Mean RT_SD IPIP6_SD PID5_SD
Total 0.220 793 601 20.34 20.48
5e66ed53de7896486f9a71f8_p8ckk 0.488 311 337 NA NA
61443e47c248a375c899a0cf_9qz3u 0.359 417 434 19.93 21.55
60d60bb9ced0566454d42750_020cx 0.349 618 227 22.37 25.55
5c43cd414fe4f800016e4983_0zp37 0.339 796 657 23.26 12.29
60d84da28f38998c15e2a4fa_xw2cl 0.333 702 376 25.12 31.38
5cbb77bf418c540001641e67_sp53j 0.324 870 547 23.39 15.09
60c5ae1d153847577f543609_u9w9l 0.302 878 950 32.58 46.24
611d69ffd22b95ca500604b2_hs0kt 0.300 855 528 24.03 15.30
6140b5b201e55201433c2383_60sey 0.299 674 303 9.79 16.05
5f7cd5bc00f2a016f22e74d2_4qe89 0.296 521 269 25.43 36.14
61564e378e974bdbe42763a2_hhufm 0.289 663 645 25.81 25.94
5c067e6329abc0000199ae45_lbvy6 0.283 579 178 14.20 17.66
60f1ce9c530c91c0a6e9dab0_1rwja 0.283 766 399 22.55 23.73
5efbc3a84e1b6416ad27a3ec_qdax2 0.281 837 2198 14.87 11.78
60e72879166656b7391b2eb2_4t9nj 0.280 897 512 19.71 25.95
6059025fd396ddcd93fe50aa_pytn9 0.279 1112 1363 20.25 23.23
613bb8d7741ae6ca04e0431b_gk94a 0.277 921 565 NA NA
5cb632bab9c4860001700fe6_dctvu 0.277 719 570 16.66 14.08
613d8489ab93bd9635c12051_zwbj3 0.272 730 269 31.39 23.54
60560efe16c645d2a9bd0daa_llvby 0.270 744 301 22.17 18.34
613c1ac80e6fc9ac43eac753_dbdal 0.266 743 396 20.36 18.70
6170adf5330be4e266555619_68ago 0.262 824 542 18.08 13.53
614b55e22ff3944a165736bb_v16g9 0.261 518 170 11.49 14.03
5eaef8702b68455d6e130595_wh7q8 0.259 707 338 14.92 12.77
606781fa7584b511b845d52c_cu5nv 0.256 550 152 NA NA
5f9aba6600cdf11f1c9b915c_at42y 0.254 788 593 22.87 20.80
5f2e1510b83f09172c62ba16_22y38 0.251 1023 1708 8.40 16.45
60feee53991033aa986fa755_lp5dv 0.251 658 330 13.83 17.35
5bdbfff9ee652a0001efca64_vvdyu 0.248 627 302 23.10 19.36
61438c83c0d9428ad7b4e4b8_t0wjb 0.248 731 336 22.27 34.83
612383021976c4e2cac79581_0fzdd 0.247 738 235 30.22 23.96
60087e9189ddb34b79b08be5_pd14q 0.245 613 202 26.64 18.85
589630d22a697d0001cfca12_buz00 0.243 850 253 18.60 30.58
5f46a8269f9d230534df642a_gql1v 0.242 1354 1254 12.23 10.61
5c437f6a4fe4f800016e3d52_a73wb 0.241 554 177 30.65 25.55
611fb1635c5f446a6bb2848d_6fcvf 0.241 990 513 37.35 28.90
5f494676e0837811e94a4007_qfzu1 0.240 940 670 16.18 21.86
5f983f3b00bbab0c4b28149a_l7hnq 0.240 646 225 11.55 21.11
5d3c6e745602310001bca8aa_6bdtb 0.238 662 495 NA NA
61645edbf8a9840feeb735b6_qot1q 0.238 1410 1660 17.48 16.72
614060a52d7c64c27ef9887c_43ppx 0.238 794 628 18.33 13.04
611e6b43e936390a65656ecd_rr1n5 0.238 578 166 27.05 19.75
5b8c8e7d0c740e00019d55c3_ejsw9 0.236 674 313 15.44 17.20
5e98b398e4a84f0eb3755e51_yy8z9 0.235 574 241 13.95 23.42
5c8fcb13dd3b360015b4dff6_dvbq9 0.234 643 327 NA NA
5f6b027e6eae971e2fa13594_5qh7m 0.233 555 180 16.80 13.21
5f304d07956a3f48639e3aa6_03qfs 0.231 698 286 26.22 27.93
615d5311a9fbf1f98aa62a5a_n2x53 0.229 495 182 18.67 20.43
610cda1332fb63830158c55d_95zqk 0.228 691 458 24.55 24.63
60ec51c51a3158a50ded8a3e_jj888 0.227 676 252 14.28 10.15
613bab4160cd2fe41f9b1a27_kepjv 0.226 587 292 20.87 22.09
60efc61b7b9f908aea94ed57_z9d93 0.224 1506 12083 26.40 17.62
612cb8236aa1cef2599e6f59_6jam2 0.221 622 751 18.93 16.50
60d1224bb86fedf1819ef67b_tvh43 0.221 790 579 17.25 26.71
6166c19fb69f46ba9749489b_6ym17 0.220 735 313 22.95 25.24
60142d79e4c5e23b4730a69e_98qsp 0.220 618 189 17.14 19.61
606f598baf5e0ffda45e9e47_qayzo 0.220 1327 1557 9.45 14.82
611ee495b12c10171def6816_9oz56 0.219 1226 928 26.97 25.54
612819d5c4f57a5fb23d03b3_sc9zu 0.218 630 467 14.06 16.49
611183c67911657721712872_sdh4z 0.218 846 529 14.36 12.28
60fa8b48cbeec7cf62317c21_glept 0.216 718 346 27.40 36.23
5c1bb460a05a64000125c522_ohxap 0.216 564 169 30.73 20.37
60feabd18109e08e594540f8_r9m9e 0.216 622 191 23.11 29.01
5e8e06cbdf10df0a811b6f29_g2tm6 0.215 610 293 NA NA
6112810f76d46087c3e2d208_51bwf 0.215 770 663 22.42 20.99
5e283bdff5bb81010fcc6cde_9zgbc 0.214 1446 1197 13.93 22.34
5d2508273854c300016c04aa_1qm3y 0.214 757 353 11.82 13.46
5b232f6838fc0c000131438c_g758c 0.212 1024 766 28.16 21.16
5e938e072a4ce9077f8157b9_egf3b 0.212 1012 781 16.75 27.39
5c42329384ef0f00011882ca_tae13 0.209 633 268 27.45 27.36
60c0fe89173c2b1f4fe2c578_a23gw 0.209 783 341 31.55 30.88
60f5fa488f7b1381d175ecb5_qwvfs 0.208 670 341 20.86 14.25
6010b797bbe6be440425665c_8qarc 0.207 816 521 11.95 15.38
605d5ffca1324379006a3599_703q5 0.207 512 139 11.48 11.99
5e6e14f53c76d23b934a67f3_ghxwv 0.206 911 2764 33.56 27.03
6139133041546096188da709_kr98d 0.206 799 465 21.71 19.13
5fcbd87e70a25c99595f012d_dkj0n 0.205 564 155 17.76 20.20
61602320c1aec7e6f8b5e3f4_adnrd 0.205 558 177 22.12 8.45
60c21f5f1389f65d2d88a3bf_lxzqo 0.204 1030 466 35.66 20.25
5f264332488e13372281c1c8_b082g 0.204 745 359 19.00 11.51
60a5404314ae1a9d5b88565d_fd4x8 0.204 763 274 15.46 15.63
61264d8a5de54eb1e42ed381_6x9p2 0.203 735 242 25.56 6.79
60009275300e1c14bf67ce83_7n71t 0.203 602 202 22.41 31.03
5fa18c4a48882016b4b143f7_kn4c0 0.202 965 984 20.97 32.36
62d6d791bb6448ae52929ebd_ncqkx 0.202 776 467 30.03 27.75
61395b705a38ede105ed7982_jkfrp 0.202 639 298 20.92 25.94
5f439ac86b921f5fac4f55ad_afpfo 0.196 884 656 22.10 16.21
60127f78468af02bcca0b7a0_lad8d 0.196 576 250 NA NA
5f97e6601f6d0e016087fc91_bac1g 0.195 649 223 20.92 19.93
614f8321333460db43b7f229_mpxhl 0.194 1645 983 NA NA
5f733312e59fd2000a74b8cf_lbgw8 0.194 972 612 12.46 16.09
5fc78656c7368206ab1c5174_4t6uk 0.192 640 306 13.37 22.43
614b23f447bd48a3801fbf95_dlqf9 0.190 739 437 16.14 18.70
6106ac34408681f3b0d07396_y2n16 0.189 613 221 27.84 24.74
60e866921fd3ba5f3007c2aa_7toga 0.189 987 1084 16.00 23.38
5ed7a7a467a98224295459ff_0l6hx 0.187 531 107 13.22 10.44
5ceef0cd45c1ec001920d548_kurpx 0.185 1062 758 22.10 22.66
610af9910c4d079d84e855f4_z8j00 0.185 764 360 25.17 16.11
5e6d381c029fa42f9f48c105_kcgfo 0.185 869 562 24.09 20.61
5f782504951dda38db3ef0b6_qsfgg 0.184 1072 798 14.89 15.58
613646a48d30fb040cf5ffed_g5q2e 0.184 702 184 14.16 15.04
614a4c9d8172751921452626_mc4l4 0.184 772 305 31.56 28.66
5fdfd04b9bf07d83b2e5f780_qcv0k 0.183 886 537 23.91 23.77
613231ee733769a9aff9c6c5_ps1tt 0.180 896 619 11.43 25.06
60d26bfdacc7efc851571e74_l8k5r 0.180 718 372 18.15 27.52
5eb3a734d249ac18a413063a_ec0gy 0.179 921 527 20.16 14.07
5f61341fe5033f0b6327fb7f_lg9vw 0.179 775 384 13.89 20.64
5ede636d21cf560bf7aeb29f_03zzd 0.178 741 532 17.33 21.62
6026f1c19366fc4778bfd055_nuccw 0.173 873 568 15.28 19.65
60ddfb3db6a71ad9ba75e387_9elmo 0.173 692 287 16.89 12.46
60d340fe43207d6a5bc89ca1_m32vm 0.173 1090 1057 14.99 18.71
5f31de301a469e11624f7948_3n2lv 0.164 783 311 19.53 26.15
611bb2a6d9678117db737236_pf2qy 0.163 904 549 12.64 12.23
6147c21c933739516362d0ee_wcaat 0.163 900 460 24.99 24.26
5ecd58b69c382e0a177e0f61_p3yev 0.162 652 279 NA NA
6067f370472020b9dc5cf71b_f2or0 0.162 861 337 16.49 17.12
61034f24da19cc56177b8b59_43rrv 0.161 805 418 28.19 13.51
614729093372b9e84bc48f8e_v4mcv 0.159 952 530 26.37 17.89
62da72f3fb2ff174670bc85f_vb5a9 0.157 786 365 18.68 13.44
5e8348450e128b0539d77ddc_z361c 0.154 820 523 9.39 24.93
5f8cbc5c355ea745e6cef2ca_5wo4f 0.152 765 406 10.78 23.51
610a9a0723dac03a63f15035_gclx4 0.151 795 312 28.52 35.59
612684c6d2cef3acadf5d77b_pv444 0.151 668 389 26.37 26.56
5d6e3252dfe6d10001dbe508_1aloo 0.150 816 379 25.32 19.75
60cefa69352cbf2549f2bf35_e13a8 0.150 688 221 24.42 25.38
5efd2964d36f63162f263795_xksu5 0.148 849 646 9.73 14.50
60ec6cb208c74c0d7e4794b9_hjay3 0.147 865 720 18.17 17.81
56cb8858edf8da000b6df354_m4821 0.146 579 211 12.57 10.43
6165e0b51a883e9db8cc7146_z0xax 0.145 816 647 23.25 14.69
5cdc88cc50f50f001a783112_tw4fr 0.118 1933 1570 23.15 13.20

Reaction Time Distribution

# RT distribution
p <- df |> 
  # filter(Date > "2022-08-12") |> 
  filter(RT < 10000) |> 
  estimate_density(select = "RT", at = c("Participant", "Block")) |>
  group_by(Participant) |>
  normalize(select = "y") |>
  ungroup() |>
  mutate(Participant = factor(Participant, levels=levels(df$Participant))) |> 
  mutate(color = ifelse(Participant %in% outliers, "red", ifelse(Participant %in% partial_outliers, "orange", "blue"))) |>
  ggplot(aes(x = x, y = y)) +
  geom_area(data = normalize(estimate_density(filter(df, RT < 10000), select = "RT"), select = "y"), alpha = 0.2) +
  geom_line(aes(color = color, group = interaction(Participant, Block), linetype = Block)) +
  geom_vline(xintercept = 125, linetype = "dashed", color = "red") +
  scale_color_manual(values = c("red" = "#F44336", "orange" = "#FF9800", "blue" = "blue"), guide = "none") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(xlim = c(0, 3000)) +
  theme_modern() +
  theme(axis.text.y = element_blank()) +
  facet_wrap(~Participant) +
  labs(y = "", x = "Reaction Time (ms)")
# p
# ggsave("figures/outliers.png", p, width=20, height=15)

# Filter out
df <- filter(df, !Participant %in% outliers)

Error Rate per Illusion Block

temp <- df |>
  group_by(Participant, Illusion_Type, Block) |>
  summarize(ErrorRate_per_block = sum(Error) / n()) |>
  ungroup() |>
  arrange(desc(ErrorRate_per_block))

temp2 <- temp |>
  filter(ErrorRate_per_block >= 0.5) |>
  group_by(Illusion_Type, Block) |>
  summarize(n = n()) |>
  arrange(desc(n), Illusion_Type) |>
  ungroup() |>
  mutate(
    n_trials = cumsum(n * 56),
    p_trials = n_trials / nrow(df)
  )

# knitr::kable(temp2)

p1 <- temp |>
  estimate_density(at = c("Illusion_Type", "Block"), method="KernSmooth") |>
  ggplot(aes(x = x, y = y)) +
  geom_line(aes(color = Illusion_Type, linetype = Block)) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(y = "Distribution", x = "Error Rate") +
  theme_modern()

p2 <- temp2 |>
  mutate(Block = fct_rev(Block)) |>
  ggplot(aes(x = Illusion_Type, y = p_trials)) +
  geom_bar(stat = "identity", aes(fill = Block)) +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(labels = scales::percent, expand = c(0, 0)) +
  labs(y = "Percentage of Trials Removed", x = "Illusion Type") +
  theme_modern() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

p1 | p2



# Drop
df <- df |>
  group_by(Participant, Illusion_Type, Block) |>
  mutate(ErrorRate_per_block = sum(Error) / n()) |>
  ungroup() |>
  filter(ErrorRate_per_block < 0.5) |>
  select(-ErrorRate_per_block)

# Drop also participant with bad second block
df <- filter(
  df,
  !(Participant %in% partial_outliers & df$Block == 2))

rm(temp, temp2)

Reaction Time per Trial

df <- df |>
  group_by(Participant, Error) |>
  mutate(Outlier = ifelse(Error == 0 & (RT < 125 | standardize(RT) > 4), TRUE, FALSE)) |>
  ungroup()

p1 <- estimate_density(df, select = "RT", at = "Participant") |>
  group_by(Participant) |>
  normalize(select = "y") |>
  ungroup() |>
  merge(df |>
    filter(Error == 0) |>
    group_by(Participant) |>
    summarize(Outlier = mean(RT) + 4 * sd(RT))) |>
  mutate(Outlier = ifelse(x >= Outlier, TRUE, FALSE)) |>
  ggplot(aes(x = x, y = y)) +
  geom_area(data = normalize(estimate_density(filter(df, RT < 10000), select = "RT"), select = "y"), alpha = 0.2) +
  geom_line(aes(color = Participant, linetype = Outlier), alpha=0.2) +
  geom_vline(xintercept = c(125), linetype = "dashed", color = "red") +
  scale_color_material_d("rainbow", guide = "none") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  guides(linetype = "none") +
  coord_cartesian(xlim = c(0, 4000)) +
  theme_modern() +
  theme(axis.text.y = element_blank()) +
  labs(y = "", x = "Reaction Time (ms)")


p2 <- df |>
  group_by(Participant) |>
  summarize(Outlier = sum(Outlier) / n()) |>
  mutate(Participant = fct_reorder(Participant, Outlier)) |>
  ggplot(aes(x = Participant, y = Outlier)) +
  geom_bar(stat = "identity", aes(fill = Participant)) +
  scale_fill_material_d("rainbow", guide = "none") +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
  see::theme_modern() +
  theme(axis.text.x = element_blank()) +
  labs(y = "Percentage of outlier trials")

p1 / p2

We removed 1244 (0.78%) outlier trials (125 ms < RT < 4 SD above mean).

df <- filter(df, Outlier == FALSE)
df$RT <- df$RT / 1000  # Convert to second for better model convergence

Participants

dfsub <- df |>
  group_by(Participant) |>
  select(Participant, Age, Sex, Education, Nationality, Ethnicity, Duration, Break_Duration, Screen_Resolution, Screen_Refresh, Device_OS, starts_with("IPIP6"), starts_with("PID5")) |>
  slice(1) |>
  ungroup()

The final sample included 128 participants (Mean age = 26.1, SD = 7.1, range: [18, 69]; Sex: 46.9% females, 52.3% males, 0.8% other; Education: Doctorate, 0.78%; Master, 12.50%; Bachelor, 46.09%; High School, 34.38%; Other, 5.47%; Prefer not to Say, 0.78%).

plot_distribution <- function(dfsub, what = "Age", title = what, subtitle = "", fill = "orange") {
  dfsub |>
    ggplot(aes_string(x = what)) +
    geom_density(fill = fill) +
    geom_vline(xintercept = mean(dfsub[[what]]), color = "red", linetype = "dashed") +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    ggtitle(title, subtitle = subtitle) +
    theme_modern() +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5),
      plot.subtitle = element_text(face = "italic", hjust = 0.5),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.y = element_blank()
    )
}

plot_waffle <- function(dfsub, what = "Nationality", title = what) {
  ggwaffle::waffle_iron(dfsub, what) |>
    # mutate(label = emojifont::fontawesome('fa-twitter')) |>
    ggplot(aes(x, y, fill = group)) +
    ggwaffle::geom_waffle() +
    # geom_point() +
    # geom_text(aes(label=label), family='fontawesome-webfont', size=4) +
    coord_equal() +
    ggtitle(title) +
    labs(fill = "") +
    theme_void() +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
p1 <- plot_distribution(dfsub, "Age", fill = "#FF9800")
p2 <- plot_distribution(dfsub, "Duration", title = "Total Duration", subtitle = "in minutes", fill = "#F44336")
p3 <- plot_distribution(dfsub, "Break_Duration", title = "Break Duration", subtitle = "in minutes", fill = "#3F51B5")

p4 <- plot_waffle(dfsub, "Sex") +
  scale_fill_manual(values = c("Male" = "#2196F3", "Female" = "#E91E63", "Other" = "#FF9800"))

p5 <- plot_waffle(dfsub, "Education") +
  scale_fill_viridis_d()

p6 <- plot_waffle(dfsub, "Nationality") +
  scale_fill_manual(values = c("South American" = "#FF5722", "Middle Eastern" = "#FF9800", "Western" = "#2196F3", "African" = "#4CAF50", "South African" = "#009688"))

p7 <- plot_waffle(dfsub, "Ethnicity") +
  scale_fill_manual(values = c("Latino" = "#FF5722", "Asian" = "#FF9800", "Caucasian" = "#2196F3", "African" = "#4CAF50", "Jewish" = "#9C27B0", "Mixed" = "#795548"))

p8 <- plot_waffle(dfsub, "Screen_Resolution", title = "Screen Resolution") +
  scale_fill_pizza_d() +
  guides(fill = "none")

p9 <- plot_waffle(dfsub, "Device_OS", title = "Device OS") +
  scale_fill_bluebrown_d()

# p10 <- plot_waffle(dfsub, "Screen_Refresh") +
#   scale_fill_viridis_d()


(p1 / p2 / p3) | (p4 / p5 / p9) | (p6 / p7 / p8)

Models

Delboeuf

Model Selection

data <- filter(df, Illusion_Type == "Delboeuf")

# mm <- model.matrix(Error ~ abs(Illusion_Strength) * Illusion_Effect, data=data)
# head(as.data.frame(mm))
best_models(data)
##                        Name   BIC R2_marginal       BF     Side  Model
## 1  DIFF_sqrt--STRENGTH_sqrt 11023      0.4205          6.62e-05 Errors
## 2  DIFF_cbrt--STRENGTH_sqrt 11026      0.4125  = 0.195 6.63e-05 Errors
## 3   DIFF_log--STRENGTH_sqrt 11038      0.4335  < 0.001 6.73e-05 Errors
## 4  DIFF_sqrt--STRENGTH_cbrt 11046      0.4345  < 0.001 6.91e-05 Errors
## 5   DIFF_sqrt--STRENGTH_log 11048      0.4116  < 0.001 6.62e-05 Errors
## 6   DIFF_sqrt--STRENGTH_log 11418      0.0710          2.70e-01     RT
## 7    DIFF_log--STRENGTH_log 11421      0.0709  = 0.175 2.66e-01     RT
## 8   DIFF_cbrt--STRENGTH_log 11422      0.0708  = 0.122 2.72e-01     RT
## 9  DIFF_sqrt--STRENGTH_sqrt 11432      0.0703  < 0.001 2.60e-01     RT
## 10      DIFF_sqrt--STRENGTH 11433      0.0702  < 0.001 2.72e-01     RT

Error Rate

Descriptive
plot_descriptive_err(data)

Model Specification
formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_delboeuf_err <- brms::brm(formula,
  data = data,
  refresh = 0
)

# parameters::parameters(model_delboeuf_err)

save(model_delboeuf_err, file="models/model_delboeuf_err.Rdata")
load("models/model_delboeuf_err.RData")
Model Visualization
p_delboeuf_err <- plot_model_err(data, model_delboeuf_err)
p_delboeuf_err

Reaction Time

Descriptive
data <- filter(df, Illusion_Type == "Delboeuf", Error == 0)

plot_descriptive_rt(data)

Model Specification
formula <- brms::bf(
  RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  sigma ~ 1 + (1 | Participant),
  beta ~ 1 + (1 | Participant),
  family = "exgaussian"
)


# brms::get_prior(formula, data = data)
# plot(seq(-5, 5, length.out=100), brms::dstudent_t(seq(-5, 5, length.out=100), df = 3, 0.1, 0.5))
priors <- c(
  # Intercept
  prior("student_t(3, 0, 0.1)", class = "Intercept"),
  prior("student_t(3, 0.35, 0.1)", class = "sd", coef="Intercept", group="Participant"),
  # Default parameters
  prior("student_t(3, 0, 3)", class = "b"),
  prior("student_t(3, 0.35, 0.5)", class = "sd", group="Participant")
  
  # prior("student_t(3, 0, 1)", class="Intercept", dpar = "sigma"),
  # prior("student_t(3, 0.1, 1)", class="sd", dpar = "sigma")
  
  # prior("student_t(3, 0.1, 0.5)", class = "ndt")
  ) |> 
  brms::validate_prior(formula, data = data)

# bayestestR::model_to_priors(model_delboeuf_rt, scale_multiply=3)

model_delboeuf_rt <- brms::brm(formula,
  data = data,
  refresh = 50,
  init=0
)

save(model_delboeuf_rt, file="models/model_delboeuf_rt.Rdata")
load("models/model_delboeuf_rt.RData")

plot_ppcheck(model_delboeuf_rt)
# performance::check_model(model_delboeuf_rt)
# parameters::parameters(model)
Model Visualization
p_delboeuf_rt <- plot_model_rt(data, model_delboeuf_rt)
p_delboeuf_rt

Visualization

p_delboeuf <- plot_all(data, p_delboeuf_err, p_delboeuf_rt)
p_delboeuf

Ebbinghaus

Model Selection

data <- filter(df, Illusion_Type == "Ebbinghaus") 
best_models(data)
##                       Name  BIC R2_marginal       BF     Side  Model
## 1  DIFF_sqrt--STRENGTH_log 8184      0.5750          8.73e-09 Errors
## 2   DIFF_log--STRENGTH_log 8185      0.5948  = 0.454 9.21e-09 Errors
## 3      DIFF_sqrt--STRENGTH 8193      0.5638  = 0.010 8.77e-09 Errors
## 4       DIFF_log--STRENGTH 8198      0.5790  = 0.001 9.33e-09 Errors
## 5  DIFF_cbrt--STRENGTH_log 8200      0.5655  < 0.001 8.83e-09 Errors
## 6      DIFF_sqrt--STRENGTH 8336      0.0943          2.31e-01     RT
## 7  DIFF_sqrt--STRENGTH_log 8338      0.0942  = 0.447 2.33e-01     RT
## 8      DIFF_cbrt--STRENGTH 8339      0.0942  = 0.320 2.32e-01     RT
## 9  DIFF_cbrt--STRENGTH_log 8340      0.0941  = 0.154 2.34e-01     RT
## 10      DIFF_log--STRENGTH 8345      0.0938  = 0.013 2.32e-01     RT

Error Rate

Descriptive
plot_descriptive_err(data)

Model Specification
formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_ebbinghaus_err <- brms::brm(formula,
  data = data,
  refresh = 0
)

# parameters::parameters(model_delboeuf_err)

save(model_ebbinghaus_err, file="models/model_ebbinghaus_err.Rdata")
load("models/model_ebbinghaus_err.RData")
Model Visualization
p_ebbinghaus_err <- plot_model_err(data, model_ebbinghaus_err)
p_ebbinghaus_err

Reaction Time

Descriptive
data <- filter(df, Illusion_Type == "Ebbinghaus", Error == 0)

plot_descriptive_rt(data)

Model Specification
formula <- brms::bf(
  RT ~ sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  sigma ~ 1 + (1 | Participant),
  beta ~ 1 + (1 | Participant),
  family = "exgaussian"
)


model_ebbinghaus_rt <- brms::brm(formula,
  data = data,
  refresh = 0,
  init = 0
)

save(model_ebbinghaus_rt, file="models/model_ebbinghaus_rt.Rdata")
load("models/model_ebbinghaus_rt.RData")

plot_ppcheck(model_ebbinghaus_rt)
# performance::check_model(model_delboeuf_rt)
# parameters::parameters(model)
Model Visualization
p_ebbinghaus_rt <- plot_model_rt(data, model_ebbinghaus_rt)
p_ebbinghaus_rt

Rod and Frame

Model Selection

data <- filter(df, Illusion_Type == "Rod-Frame")
best_models(data)
##                        Name   BIC R2_marginal       BF     Side  Model
## 1    DIFF_log--STRENGTH_log 10354      0.4875          0.000146 Errors
## 2   DIFF_sqrt--STRENGTH_log 10357      0.4893  = 0.208 0.000148 Errors
## 3   DIFF_log--STRENGTH_sqrt 10367      0.4845  = 0.002 0.000146 Errors
## 4  DIFF_sqrt--STRENGTH_sqrt 10370      0.4865  < 0.001 0.000147 Errors
## 5   DIFF_log--STRENGTH_cbrt 10374      0.4905  < 0.001 0.000149 Errors
## 6    DIFF_log--STRENGTH_log  8688      0.0900          0.076184     RT
## 7   DIFF_log--STRENGTH_sqrt  8689      0.0899  = 0.435 0.076600     RT
## 8        DIFF_log--STRENGTH  8695      0.0895  = 0.028 0.079896     RT
## 9   DIFF_log--STRENGTH_cbrt  8700      0.0891  = 0.002 0.077029     RT
## 10  DIFF_sqrt--STRENGTH_log  8701      0.0891  = 0.001 0.077129     RT

Error Rate

Descriptive
plot_descriptive_err(data)

Model Specification
formula <- brms::bf(
  Error ~ logmod(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + logmod(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_rodframe_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 1585.6 seconds.
## Chain 1 finished in 1610.1 seconds.
## Chain 2 finished in 1616.1 seconds.
## Chain 4 finished in 1618.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1607.6 seconds.
## Total execution time: 1618.8 seconds.

# parameters::parameters(model_delboeuf_err)

save(model_rodframe_err, file="models/model_rodframe_err.Rdata")
load("models/model_rodframe_err.RData")
Model Visualization
p_rodframe_err <- plot_model_err(data, model_rodframe_err)
p_rodframe_err

Reaction Time

Descriptive
data <- filter(df, Illusion_Type == "Rod-Frame", Error == 0)

plot_descriptive_rt(data)

Model Specification
formula <- brms::bf(
  RT ~ logmod(Illusion_Difference) * logmod(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  sigma ~ 1 + (1 | Participant),
  beta ~ 1 + (1 | Participant),
  family = "exgaussian"
)


model_rodframe_rt <- brms::brm(formula,
  data = data,
  refresh = 0,
  init = 0
)

save(model_rodframe_rt, file="models/model_rodframe_rt.Rdata")
load("models/model_rodframe_rt.RData")

plot_ppcheck(model_rodframe_rt)
# performance::check_model(model_delboeuf_rt)
# parameters::parameters(model)
Model Visualization
p_rodframe_rt <- plot_model_rt(data, model_rodframe_rt)
p_rodframe_rt

Vertical-Horizontal

Model Selection

data <- filter(df, Illusion_Type == "Vertical-Horizontal") 
best_models(data)
##                        Name   BIC R2_marginal       BF  Side  Model
## 1  DIFF_sqrt--STRENGTH_sqrt 10047       0.579          0.128 Errors
## 2  DIFF_cbrt--STRENGTH_sqrt 10054       0.578  = 0.027 0.129 Errors
## 3   DIFF_log--STRENGTH_sqrt 10061       0.581  < 0.001 0.129 Errors
## 4       DIFF--STRENGTH_sqrt 10071       0.582  < 0.001 0.129 Errors
## 5  DIFF_sqrt--STRENGTH_cbrt 10143       0.586  < 0.001 0.129 Errors
## 6   DIFF_log--STRENGTH_sqrt  7234       0.128          0.875     RT
## 7  DIFF_sqrt--STRENGTH_sqrt  7234       0.128  = 0.956 0.888     RT
## 8       DIFF--STRENGTH_sqrt  7237       0.128  = 0.189 0.871     RT
## 9  DIFF_cbrt--STRENGTH_sqrt  7238       0.128  = 0.089 0.893     RT
## 10  DIFF_log--STRENGTH_cbrt  7247       0.127  = 0.001 0.880     RT

Error Rate

Descriptive
plot_descriptive_err(data)

Model Specification
formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_verticalhorizontal_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 1850.4 seconds.
## Chain 4 finished in 1869.3 seconds.
## Chain 2 finished in 1916.9 seconds.
## Chain 3 finished in 1929.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1891.5 seconds.
## Total execution time: 1929.7 seconds.

# parameters::parameters(model_verticalhorizontal_err)

save(model_verticalhorizontal_err, file="models/model_verticalhorizontal_err.Rdata")
load("models/model_verticalhorizontal_err.RData")
Model Visualization
p_verticalhorizontal_err <- plot_model_err(data, model_verticalhorizontal_err)
p_verticalhorizontal_err

Reaction Time

Descriptive
data <- filter(df, Illusion_Type == "Vertical-Horizontal", Error == 0)

plot_descriptive_rt(data)

Model Specification
formula <- brms::bf(
  RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ t2(Illusion_Difference * Illusion_Strength) +
  #   (1 | Participant),
  family = "exgaussian"
)


model_verticalhorizontal_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_verticalhorizontal_rt, file="models/model_verticalhorizontal_rt.Rdata")
load("models/model_verticalhorizontal_rt.RData")

plot_ppcheck(model_verticalhorizontal_rt)
# performance::check_model(model_verticalhorizontal_rt)
# parameters::parameters(model)
Model Visualization
p_verticalhorizontal_rt <- plot_model_rt(data, model_verticalhorizontal_rt)
p_verticalhorizontal_rt

Zöllner

Model Selection

data <- filter(df, Illusion_Type == "Zöllner")
best_models(data)
##                        Name   BIC R2_marginal       BF    Side  Model
## 1       DIFF_cbrt--STRENGTH  9367      0.4069          0.71416 Errors
## 2        DIFF_log--STRENGTH  9375      0.4154  = 0.013 0.71320 Errors
## 3       DIFF_sqrt--STRENGTH  9376      0.4170  = 0.009 0.71262 Errors
## 4            DIFF--STRENGTH  9521      0.4468  < 0.001 0.70663 Errors
## 5  DIFF_cbrt--STRENGTH_sqrt  9772      0.3878  < 0.001 0.70402 Errors
## 6       DIFF_cbrt--STRENGTH 13383      0.0690          0.00455     RT
## 7        DIFF_log--STRENGTH 13384      0.0689  = 0.545 0.00467     RT
## 8       DIFF_sqrt--STRENGTH 13396      0.0683  = 0.001 0.00483     RT
## 9  DIFF_cbrt--STRENGTH_sqrt 13420      0.0670  < 0.001 0.00442     RT
## 10  DIFF_log--STRENGTH_sqrt 13420      0.0669  < 0.001 0.00455     RT

Error Rate

Descriptive
plot_descriptive_err(data)

Model Specification
formula <- brms::bf(
  Error ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
    (1 + cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_zollner_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 1926.9 seconds.
## Chain 1 finished in 1964.4 seconds.
## Chain 2 finished in 1973.6 seconds.
## Chain 4 finished in 1980.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1961.3 seconds.
## Total execution time: 1980.3 seconds.

# parameters::parameters(model_zollner_err)

save(model_zollner_err, file="models/model_zollner_err.Rdata")
load("models/model_zollner_err.RData")
Model Visualization
p_zollner_err <- plot_model_err(data, model_zollner_err)
p_zollner_err

Reaction Time

Descriptive
data <- filter(df, Illusion_Type == "Zöllner", Error == 0)

plot_descriptive_rt(data)

Model Specification
formula <- brms::bf(
  RT ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ t2(Illusion_Difference * Illusion_Strength) +
  #   (1 | Participant),
  family = "exgaussian"
)


model_zollner_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_zollner_rt, file="models/model_zollner_rt.Rdata")
load("models/model_zollner_rt.RData")

plot_ppcheck(model_zollner_rt)
# performance::check_model(model_zollner_rt)
# parameters::parameters(model)
Model Visualization
p_zollner_rt <- plot_model_rt(data, model_zollner_rt)
p_zollner_rt

White

data <- filter(df, Illusion_Type == "White")
best_models(data)
##                        Name  BIC R2_marginal       BF    Side  Model
## 1  DIFF_sqrt--STRENGTH_sqrt 5821      0.7677          0.00175 Errors
## 2       DIFF--STRENGTH_sqrt 5827      0.7809  = 0.045 0.00165 Errors
## 3  DIFF_cbrt--STRENGTH_sqrt 5834      0.7644  = 0.002 0.00182 Errors
## 4   DIFF_log--STRENGTH_sqrt 5861      0.7615  < 0.001 0.00193 Errors
## 5   DIFF_sqrt--STRENGTH_log 5862      0.7970  < 0.001 0.00200 Errors
## 6        DIFF_log--STRENGTH 9590      0.0903          0.01015     RT
## 7       DIFF_cbrt--STRENGTH 9600      0.0897  = 0.008 0.01027     RT
## 8       DIFF_sqrt--STRENGTH 9610      0.0892  < 0.001 0.01035     RT
## 9   DIFF_log--STRENGTH_sqrt 9644      0.0866  < 0.001 0.01129     RT
## 10           DIFF--STRENGTH 9648      0.0869  < 0.001 0.01048     RT

Error Rate

Descriptive
plot_descriptive_err(data)

Model Specification
formula <- brms::bf(
  Error ~ Illusion_Difference * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + Illusion_Difference * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_white_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Start sampling
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 2677.9 seconds.
## Chain 4 finished in 2991.8 seconds.
## Chain 3 finished in 3058.0 seconds.
## Chain 2 finished in 3117.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 2961.3 seconds.
## Total execution time: 3117.8 seconds.

# parameters::parameters(model_white_err)

save(model_white_err, file="models/model_white_err.Rdata")
load("models/model_white_err.RData")
Model Visualization
p_white_err <- plot_model_err(data, model_white_err)
p_white_err

Reaction Time

Descriptive
data <- filter(df, Illusion_Type == "White", Error == 0)

plot_descriptive_rt(data)

Model Specification
formula <- brms::bf(
  RT ~ Illusion_Difference * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ t2(Illusion_Difference * Illusion_Strength) +
  #   (1 | Participant),
  family = "exgaussian"
)


model_white_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_white_rt, file="models/model_white_rt.Rdata")
load("models/model_white_rt.RData")

plot_ppcheck(model_white_rt)
# performance::check_model(model_white_rt)
# parameters::parameters(model)
Model Visualization
p_white_rt <- plot_model_rt(data, model_white_rt)
p_white_rt

Müller-Lyer

Model Selection

data <- filter(df, Illusion_Type == "Müller-Lyer")
best_models(data)
##                        Name  BIC R2_marginal       BF     Side  Model
## 1  DIFF_sqrt--STRENGTH_sqrt 7273       0.732          5.25e-17 Errors
## 2   DIFF_log--STRENGTH_sqrt 7273       0.730  = 0.738 4.28e-17 Errors
## 3       DIFF--STRENGTH_sqrt 7291       0.730  < 0.001 4.32e-17 Errors
## 4  DIFF_cbrt--STRENGTH_sqrt 7292       0.734  < 0.001 6.54e-17 Errors
## 5    DIFF_log--STRENGTH_log 7299       0.752  < 0.001 5.76e-17 Errors
## 6       DIFF_sqrt--STRENGTH 9593       0.154          1.64e-09     RT
## 7       DIFF_cbrt--STRENGTH 9593       0.154  = 0.953 1.62e-09     RT
## 8  DIFF_cbrt--STRENGTH_sqrt 9602       0.153  = 0.007 6.02e-10     RT
## 9  DIFF_sqrt--STRENGTH_sqrt 9604       0.153  = 0.004 6.19e-10     RT
## 10       DIFF_log--STRENGTH 9605       0.153  = 0.002 1.86e-09     RT

Error Rate

Descriptive
plot_descriptive_err(data, side = "updown")

Model Specification
formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_mullerlyer_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 2054.6 seconds.
## Chain 1 finished in 2112.0 seconds.
## Chain 4 finished in 2506.1 seconds.
## Chain 2 finished in 2610.9 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 2320.9 seconds.
## Total execution time: 2611.3 seconds.

# parameters::parameters(model_white_err)

save(model_mullerlyer_err, file="models/model_mullerlyer_err.Rdata")
load("models/model_mullerlyer_err.RData")
Model Visualization
p_mullerlyer_err <- plot_model_err(data, model_mullerlyer_err)
p_mullerlyer_err

Reaction Time

Descriptive
data <- filter(df, Illusion_Type == "Müller-Lyer", Error == 0)

plot_descriptive_rt(data, side = "updown")

Model Specification
formula <- brms::bf(
  RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ t2(Illusion_Difference * Illusion_Strength) +
  #   (1 | Participant),
  family = "exgaussian"
)


model_mullerlyer_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_mullerlyer_rt, file="models/model_mullerlyer_rt.Rdata")
# parameters::parameters(model_white_err)
load("models/model_mullerlyer_rt.RData")

plot_ppcheck(model_mullerlyer_rt)
# performance::check_model(model_mullerlyer_rt)
# parameters::parameters(model)
Model Visualization
p_mullerlyer_rt <- plot_model_rt(data, model_mullerlyer_rt)
p_mullerlyer_rt

Ponzo

Model Selection

data <- filter(df, Illusion_Type == "Ponzo")
best_models(data)
##                        Name   BIC R2_marginal       BF     Side  Model
## 1       DIFF_cbrt--STRENGTH  8250       0.563          0.735137 Errors
## 2       DIFF_sqrt--STRENGTH  8254       0.570  = 0.135 0.737259 Errors
## 3        DIFF_log--STRENGTH  8304       0.590  < 0.001 0.742920 Errors
## 4            DIFF--STRENGTH  8339       0.599  < 0.001 0.745206 Errors
## 5  DIFF_cbrt--STRENGTH_sqrt  8368       0.586  < 0.001 0.761796 Errors
## 6       DIFF_cbrt--STRENGTH 10267       0.118          0.000664     RT
## 7       DIFF_sqrt--STRENGTH 10279       0.117  = 0.003 0.000631     RT
## 8        DIFF_log--STRENGTH 10325       0.115  < 0.001 0.000579     RT
## 9            DIFF--STRENGTH 10357       0.114  < 0.001 0.000566     RT
## 10 DIFF_cbrt--STRENGTH_sqrt 10362       0.113  < 0.001 0.001393     RT

Error Rate

Descriptive
plot_descriptive_err(data, side = "updown")

Model Specification
formula <- brms::bf(
  Error ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
    (1 + cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_ponzo_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 2044.8 seconds.
## Chain 4 finished in 2077.4 seconds.
## Chain 3 finished in 2093.1 seconds.
## Chain 1 finished in 2096.9 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 2078.1 seconds.
## Total execution time: 2097.2 seconds.

# parameters::parameters(model_ponzo_err)

save(model_ponzo_err, file="models/model_ponzo_err.Rdata")
load("models/model_ponzo_err.RData")
Model Visualization
p_ponzo_err <- plot_model_err(data, model_ponzo_err)
p_ponzo_err

Reaction Time

Descriptive
data <- filter(df, Illusion_Type == "Ponzo", Error == 0)

plot_descriptive_rt(data, side = "updown")

Model Specification
formula <- brms::bf(
  RT ~ cbrtmod(Illusion_Difference) * abs(Illusion_Strength) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ t2(Illusion_Difference * Illusion_Strength) +
  #   (1 | Participant),
  family = "exgaussian"
)


model_ponzo_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_ponzo_rt, file="models/model_ponzo_rt.Rdata")
# parameters::parameters(model_white_err)
load("models/model_ponzo_rt.RData")

plot_ppcheck(model_ponzo_rt)
# performance::check_model(model_ponzo_rt)
# parameters::parameters(model)
Model Visualization
p_ponzo_rt <- plot_model_rt(data, model_ponzo_rt)
p_ponzo_rt

Poggendorff

Model Selection

data <- filter(df, Illusion_Type == "Poggendorff")
best_models(data)
##                        Name  BIC R2_marginal       BF   Side  Model
## 1  DIFF_cbrt--STRENGTH_sqrt 8259      0.5190          0.6545 Errors
## 2  DIFF_sqrt--STRENGTH_sqrt 8281      0.5193  < 0.001 0.6607 Errors
## 3       DIFF_cbrt--STRENGTH 8325      0.4986  < 0.001 0.6401 Errors
## 4       DIFF_sqrt--STRENGTH 8348      0.4973  < 0.001 0.6463 Errors
## 5   DIFF_log--STRENGTH_sqrt 8380      0.5184  < 0.001 0.6746 Errors
## 6       DIFF_cbrt--STRENGTH 8461      0.0778          0.0251     RT
## 7       DIFF_sqrt--STRENGTH 8480      0.0766  < 0.001 0.0246     RT
## 8  DIFF_cbrt--STRENGTH_sqrt 8520      0.0740  < 0.001 0.0268     RT
## 9        DIFF_log--STRENGTH 8534      0.0734  < 0.001 0.0233     RT
## 10 DIFF_sqrt--STRENGTH_sqrt 8540      0.0728  < 0.001 0.0263     RT

Error Rate

Descriptive
plot_descriptive_err(data, side = "updown")

Model Specification
formula <- brms::bf(
  Error ~ cbrtmod(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + cbrtmod(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_poggendorff_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 3086.8 seconds.
## Chain 3 finished in 3150.5 seconds.
## Chain 2 finished in 3172.7 seconds.
## Chain 1 finished in 3187.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 3149.4 seconds.
## Total execution time: 3187.6 seconds.

# parameters::parameters(model_poggendorff_err)

save(model_poggendorff_err, file="models/model_poggendorff_err.Rdata")
load("models/model_poggendorff_err.RData")
Model Visualization
p_poggendorff_err <- plot_model_err(data, model_poggendorff_err)
p_poggendorff_err

Reaction Time

Descriptive
data <- filter(df, Illusion_Type == "Poggendorff", Error == 0)

plot_descriptive_rt(data, side = "updown")

Model Specification
formula <- brms::bf(
  RT ~ cbrtmod(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ 1 + (1 | Participant),
  family = "exgaussian"
)


model_poggendorff_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_ponzo_rt, file="models/model_poggendorff_rt.Rdata")
load("models/model_poggendorff_rt.RData")

plot_ppcheck(model_poggendorff_rt)
# performance::check_model(model_poggendorff_rt)
# parameters::parameters(model)
Model Visualization
p_poggendorff_rt <- plot_model_rt(data, model_poggendorff_rt)
p_poggendorff_rt

Contrast

Model Selection

data <- filter(df, Illusion_Type == "Contrast")
best_models(data)
##                        Name   BIC R2_marginal       BF     Side  Model
## 1  DIFF_sqrt--STRENGTH_sqrt  7837       0.610          4.70e-19 Errors
## 2       DIFF--STRENGTH_sqrt  7842       0.621  = 0.084 4.55e-19 Errors
## 3  DIFF_cbrt--STRENGTH_sqrt  7847       0.607  = 0.008 5.12e-19 Errors
## 4   DIFF_sqrt--STRENGTH_log  7859       0.658  < 0.001 8.41e-19 Errors
## 5        DIFF--STRENGTH_log  7864       0.675  < 0.001 8.08e-19 Errors
## 6  DIFF_sqrt--STRENGTH_sqrt 10454       0.122          2.42e-02     RT
## 7  DIFF_cbrt--STRENGTH_sqrt 10454       0.122  = 0.824 2.43e-02     RT
## 8   DIFF_log--STRENGTH_sqrt 10457       0.122  = 0.211 2.45e-02     RT
## 9       DIFF--STRENGTH_sqrt 10462       0.122  = 0.013 2.44e-02     RT
## 10      DIFF_sqrt--STRENGTH 10479       0.121  < 0.001 1.94e-02     RT

Error Rate

Descriptive
data <- filter(df, Illusion_Type == "Contrast")

plot_descriptive_err(data, side = "updown")

Model Specification
formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 + sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_contrast_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 3483.4 seconds.
## Chain 2 finished in 3506.1 seconds.
## Chain 4 finished in 3551.5 seconds.
## Chain 3 finished in 3555.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 3524.0 seconds.
## Total execution time: 3555.6 seconds.

# parameters::parameters(model_contrast_err)

save(model_contrast_err, file="models/model_contrast_err.Rdata")
load("models/model_contrast_err.RData")
Model Visualization
p_contrast_err <- plot_model_err(data, model_contrast_err)
p_contrast_err

Reaction Time

Descriptive
data <- filter(df, Illusion_Type == "Contrast", Error == 0)

plot_descriptive_rt(data, side = "updown")

Model Specification
formula <- brms::bf(
  RT ~ sqrt(Illusion_Difference) * sqrt(abs(Illusion_Strength)) * Illusion_Effect +
    (1 | Participant),
  # sigma ~ 1 + (1 | Participant),
  family = "exgaussian"
)


model_contrast_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)

save(model_contrast_rt, file="models/model_contrast_rt.Rdata")
load("models/model_contrast_rt.RData")

plot_ppcheck(model_contrast_rt)
# performance::check_model(model_contrast_rt)
# parameters::parameters(model)
Model Visualization
p_contrast_rt <- plot_model_rt(data, model_contrast_rt)
p_contrast_rt

Individual Scores

Reliability

random <- rbind(
  extract_random(model_delboeuf_err, "Delboeuf"),
  # extract_random(model_delboeuf_rt, "Delboeuf")
  extract_random(model_ebbinghaus_err, "Ebbinghaus"),
  # extract_random(model_ebbinghaus_rt, "Ebbinghaus"),
  extract_random(model_rodframe_err, "Rod-Frame"),
  # extract_random(model_rodframe_rt, "Rod-Frame"),
  extract_random(model_verticalhorizontal_err, "Vertical-Horizontal"),
  # extract_random(model_verticalhorizontal_rt, "Vertical-Horizontal"),
  extract_random(model_zollner_err, "Zöllner"),
  # extract_random(model_zollner_rt, "Zöllner"),
  extract_random(model_white_err, "White"),
  # extract_random(model_white_rt, "White"),
  extract_random(model_mullerlyer_err, "Müller-Lyer"),
  # extract_random(model_mullerlyer_rt, "Müller-Lyer"),
  extract_random(model_ponzo_err, "Ponzo"),
  # extract_random(model_ponzo_rt, "Ponzo"),
  extract_random(model_poggendorff_err, "Poggendorff"),
  # extract_random(model_poggendorff_rt, "Poggendorff"),
  extract_random(model_contrast_err, "Contrast")
  # extract_random(model_contrast_rt, "Contrast")
) |> 
  filter(
    # !str_detect(Parameter, "Intercept_"),
    !str_detect(Parameter, "Cong_"),
    !str_detect(Parameter, "Null_"))


random |>
  # group_by(Illusion_Type, Parameter) |> 
  # standardize(select="Value") |> 
  # ungroup() |> 
  # mutate(order = mean(Value)) |> 
  # ungroup() |> 
  # mutate(Participant = fct_reorder(Participant, order)) |> 
  ggplot(aes(y = Participant, x = Value)) +
  # ggdist::stat_slabinterval(aes(fill = Illusion_Type, color=Illusion_Type), alpha=0.5) +
  ggdist::stat_pointinterval(aes(fill = Illusion_Type, color=Illusion_Type), alpha=0.5) +
  # coord_cartesian(xlim=c(-4, 4)) +
  facet_wrap(~Parameter, scales = "free_x") +
  theme(axis.text.y = element_blank())



random <- random |>
  mutate(Parameter = paste0(Illusion_Type, "_", Parameter),
         Parameter = format_illusionName(Parameter)) |>
  group_by(Participant, Parameter) |>
  summarize(Score = median(Value), .groups = "drop") |>
  ungroup() |>
  group_by(Parameter) |> 
  mutate(Score = Score / sd(Score)) |>  # Scale
  ungroup()  |> 
  pivot_wider(names_from = Parameter, values_from = Score)

# plot(estimate_density(select(random, -Participant)))

Relationship with empirical descriptors

# Average when incongruent
empirical <- df |>
  filter(Illusion_Effect == "Incongruent") |> 
  group_by(Participant, Illusion_Type) |>
  summarize(
    # n = n(),
    Error_General = sum(Error, na.rm=TRUE) / n(),
    RT_Mean_General = mean(RT, na.rm=TRUE),
    RT_SD_General = sd(RT, na.rm=TRUE),
  ) |>
  ungroup() |> 
  mutate(Illusion_Type = format_illusionName(Illusion_Type))  |> 
  pivot_wider(names_from = c("Illusion_Type"), values_from = c("Error_General", "RT_Mean_General", "RT_SD_General")) |> 
  arrange(Participant)

# Average for 3 higher level of illusion strength
empirical <- df |>
  filter(Illusion_Effect == "Incongruent") |> 
  group_by(Participant, Illusion_Type, Illusion_Strength) |>
  summarize(
    # n = n(),
    Error_Strong = sum(Error, na.rm=TRUE) / n(),
    RT_Mean_Strong = mean(RT, na.rm=TRUE),
    RT_SD_Strong = sd(RT, na.rm=TRUE),
  ) |>
  slice((n()-2):n()) |> 
  group_by(Participant, Illusion_Type) |>
  summarize(
    # n = n(),
    Error_Strong = mean(Error_Strong, na.rm=TRUE),
    RT_Mean_Strong = mean(RT_Mean_Strong, na.rm=TRUE),
    RT_SD_Strong = mean(RT_SD_Strong, na.rm=TRUE),
  ) |> 
  ungroup() |> 
  mutate(Illusion_Type = format_illusionName(Illusion_Type)) |> 
  pivot_wider(names_from = c("Illusion_Type"), values_from = c("Error_Strong", "RT_Mean_Strong", "RT_SD_Strong")) |> 
  arrange(Participant) |> 
  cbind(empirical)

# Correlation
r <- data.frame()
for(ill in c("Delboeuf", "Ebbinghaus", 
             # "RodFrame", 
             "VerticalHorizontal", 
             # "Zollner", 
             # "White", 
             "MullerLyer", 
             "Ponzo", 
             # "Poggendorff", 
             "Contrast")) {
  r <- correlation::correlation(
    data = random |> 
      arrange(Participant) |> 
      data_select(select=ill, regex=TRUE),
    data2 = empirical |>
      data_select(select=ill, regex=TRUE) 
    ) |> 
    mutate(Illusion_Type = ill,
           Parameter1 = clean_name(Parameter1),
           Parameter2 = clean_name(Parameter2)) |> 
    select(Illusion_Type, Parameter1, Parameter2, r, p) |> 
    rbind(r)
}

r |>
  ggplot(aes(x = Parameter2, y = Parameter1)) +
  geom_tile(aes(fill = r, alpha = p)) +
  scale_alpha_continuous(range = c(1, 0)) +
  scale_fill_gradient2(low = "#2196F3", mid = "white", high = "#F44336", midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Correlation", guide = "legend") +
  facet_wrap(~Illusion_Type, scales = "free")

Structure

Correlation

cor <- correlation::correlation(random, redundant = TRUE)

p_data <- cor |>
  mutate(
    Parameter1 = clean_name(Parameter1),
    Parameter2 = clean_name(Parameter2)
  ) |>
  correlation::cor_sort(hclust_method = "ward.D2") |>
  cor_lower() |>
  mutate(
    Text = insight::format_value(r, zap_small = TRUE, digits = 3),
    Text = str_replace(str_remove(Text, "^0+"), "^-0+", "-"),
    Parameter2 = fct_rev(Parameter2)
  )


p_data |>
  ggplot(aes(x = Parameter2, y = Parameter1)) +
  geom_tile(aes(fill = r)) +
  # geom_text(aes(label = Text), size = 2) +
  scale_fill_gradient2(low = "#2196F3", mid = "white", high = "#F44336", midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Correlation", guide = "legend") +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_discrete(expand = c(0, 0)) +
  labs(title = "Correlation Matrix", x = NULL, y = NULL) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

Factor Analysis

First Level

data <- select(random, -Participant)

rez <- parameters::n_factors(data)
plot(rez)


fa <- parameters::factor_analysis(data,
  cor = as.matrix(cor),
  n = 8,
  rotation = "oblimin",
  fm = "mle",
  sort = TRUE
)

# psych::omega(data, nfactors=3)
# fa <- psych::fa.multi(as.matrix(cor), nfactors = 3, nfact2 = 1, n.obs = nrow(random))
# psych::fa.multi.diagram(fa)

Second Level

data <- predict(fa)

rez <- parameters::n_factors(data, rotation = "oblimin")
plot(rez)

fa2 <- parameters::factor_analysis(data,
  n = 2,
  rotation = "varimax",
  fm = "mle",
  sort = TRUE
)

# psych::schmid(attributes(fa)$model, nfactors = 1)

Visualization

library(ggraph)

dat <- rbind(
  attributes(fa)$loadings_long,
  attributes(fa2)$loadings_long |>
    mutate(Component = str_replace(Component, "ML", "G"))
) |>
  tidygraph::as_tbl_graph() |>
  tidygraph::activate("nodes")


# Get order of indices
# idx <- correlation::cor_sort(cor, hclust_method = "ward.D2") |>
#   summary() |>
#   pull(Parameter)
# idx <- sort(fa$Variable)
idx <- fa$Variable

get_x <- function(name, x) {
  case_when(
    name == "G1" ~ 0/1,
    name == "G2" ~ 1/1,
    name == "ML1" ~ 0/4,
    name == "ML2" ~ 1/4,
    name == "ML3" ~ 2/4,
    name == "ML4" ~ 3/4,
    name == "ML5" ~ 4/4,
    TRUE ~ x)
}

dat |>
  mutate(
    x = sapply(as.data.frame(dat)$name, function(x) {
      if (x %in% idx) {
        (which(idx == x) - 1) / (length(idx) - 1)
      } else {
        NA
      }
    })
  ) |>
  mutate(
    x = get_x(name, x),
    y = case_when(
    str_detect(name, "G") ~ 2,
    str_detect(name, "ML") ~ 1,
    TRUE ~ 0),
    angle = ifelse(name %in% idx, 90, 0),
    hjust = ifelse(name %in% idx, 1.01, 0),
    name = clean_name(name),
    Illusion_Type = str_split_fixed(name, " - ", n=2)[, 1]) |> 
  tidygraph::activate("edges") |> 
  filter(abs(Loading) > 0.2) |> 
  ggraph(layout = "nicely") +
  geom_edge_link(aes(edge_width = Loading, edge_colour = as.factor(sign(Loading)))) +
  geom_node_point(aes(color = Illusion_Type), size=3) +
  geom_node_text(aes(label = name, angle = angle, hjust=hjust)) +
  scale_y_continuous(expand = expansion(c(.6, .1))) +
  scale_edge_width_continuous(range=c(0.1, 1)) +
  scale_edge_color_manual(values=c("#F44336", "#2196F3")) +
  scale_color_manual(values=c("Delboeuf" = "#F44336", Ebbinghaus = "#E91E63",
                              "Rod-Frame" = "#9C27B0")) +
  guides(edge_colour = "none") +
  ggraph::theme_graph()

SEM

library(lavaan)
library(tidySEM)


names <- names(select(random, -Participant))
data <- setNames(as.data.frame(str_split_fixed(names, "_", 3)), c("Illusion_Type", "Parameter", "Model"))
data$Name <- names

# Model 1
lvl1 <- paste(data$Illusion_Type, "=~", data$Name)
lvl2 <- paste("G =~", unique(data$Illusion_Type))
resid <- as.data.frame(t(utils::combn(x=unique(data$Illusion_Type), 2)))
resid <- paste(resid$V1, "~~", resid$V2)[1:2]



m1 <- sem(model = paste0(c(lvl1, lvl2, resid), collapse = "\n"),
           data  = random)


# Model 2
lvl1 <- paste(data$Model, "=~", data$Name)
lvl2 <- paste("G =~", unique(data$Model))

m2 <- sem(model = paste0(c(lvl1, lvl2), collapse = "\n"),
           data  = random)

# Model 3
lvl1 <- paste(data$Parameter, "=~", data$Name)
lvl2 <- paste("G =~", unique(data$Parameter))

m3 <- sem(model = paste0(c(lvl1, lvl2), collapse = "\n"),
           data  = random)

anova(m1, m2, m3)

summary(m1, standardize = TRUE)

tidySEM::graph_sem(model = m1)

# as_tbl_graph(fit, standardize = TRUE)
# Save
dfsub <- dfsub |> 
  merge(random, by = "Participant") |> 
  cbind(data, predict(fa2, names=c("G1", "G2")))

write.csv("../data/study3.csv")